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ABSTRACT 

We use a self-consistent Monte Carlo treatment of stellar dynamics to investigate 
black hole binaries that are dynamically ejected from globular clusters to determine 
if they will be gravitational wave sources. We find that many of the ejected binaries 
have initially short periods and will merge within a Hubble time due to gravitational 
wave radiation. Thus they are potential sources for ground-based gravitational wave 
detectors. We estimate the yearly detection rate for current and advanced ground- 
based detectors and find a modest enhancement over the rate predicted for binaries 
produced by pure stellar evolution in galactic fields. We also find that many of the 
ejected binaries will pass through the longer wavelength Laser Interferometer Space 
Antenna (LISA) band and may be individually resolvable. We find a low probability 
that the Galaxy will contain a binary in the LISA band during its three-year mis- 
sion. Some such binaries may, however, be detectable at Mpc distances implying that 
there may be resolvable stellar-mass LISA sources beyond our Galaxy. We conclude 
that globular clusters have a significant effect on the detection rate of ground-based 
detectors and may produce interesting LISA sources in local group galaxies. 

Key words: galaxies: star clusters - gravitational waves - stellar dynamics - binaries: 
close - globular clusters: general - methods: N-body 



1 INTRODUCTION 

Gravitational waves, metric perturbations in space-time 
caused by time-varying mass-energy distributions, offer a 
new window on the universe that is independent of the elec- 
tromagnetic spectrum. In particular the inspiral and merger 
of stellar-mass compact binaries (binaries consisting of neu- 
tron stars (NSs) and/or black holes (BHs)) are predicted 
to be major burst sources for the high-frequency ground- 
based gravitational wave detectors Virgo and LIGO. The 
current generation of these detectors should be able to detect 
such events out the the Virgo cluster (Abbott et al. 2005, 
2006, 2010) while upgrades to these detectors (e.g. advanced 
LIGO) should lead to detections at cosmologically signifi- 
cant distances. At larger separations compact binaries, in- 
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eluding white dwarf (WD) binaries, within the Galaxy may 
be detected by the low-frequency space-based gravitational 
wave detector LISA planned for 2018-2020 (Hils, Bender 
& Webbink 1990; Benacquista 2001a; Nelemans, Yungelson 
& Portegies Zwart 2001; Belczynski, Benacquista & Bulik 
2008) . WD- WD binaries are expected to be plentiful and will 
manifest as a confusion limited noise source (Evans, Iben & 
Smarr 1987; Hils, Bender & Webbink 1990; Nelemans, Yun- 
gelson & Portegies Zwart 2001; Timpano, Rubbo & Cornish 
2006; Ruiter et al. 2010) whereas the rarer NS-NS, NS-BH 
and BH-BH binaries may be individually resolved. 

Many of the sources producing gravitational waves will 
not emit significant amounts of electromagnetic radiation 
(particularly NS-NS and BH-BH binaries) and it is necessary 
to perform population synthesis models in order to constrain 
event rates and produce templates for gravitational wave 
detectors. Significant efforts have been made to constrain 
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the population of compact binaries in galactic fields (e.g. 
Portegies Zwart & Yungelson 1998; Fryer, Woosley & Hart- 
mann 1999; Belczynski, Kalogera & Bulik 2002; Belczynski 
et al. 2007) and these studies have found that detection rates 
should be dominated by NS-NS inspirals. Although BH-BH 
binaries are more massive and can be detected with a higher 
signal-to-noise ratio at larger distances, NSs are both more 
plentiful due to the power-law shape of the initial mass func- 
tion (IMF) and the progenitors of NS-NS binaries are less 
likely to undergo mass-transfer leading to a merger than are 
their BH-BH counterparts. In particular Belczynski et al. 
(2007) determined advanced LIGO detection rates of ~ 20 
yr~^ for NS-NS binaries, ~ 2 yr^^ for BH-BH binaries and 
only ~ 1 yr~^ for NS-BH binaries. Belczynski, Benacquista 
& Bulik (2008) have similarly concluded that the number 
of resolved stellar-mass detections by LISA, although small, 
will also be dominated by NS-NS binaries. 

Although BH-BH binaries may be rare in galactic fields, 
it is possible to produce such binaries through dynamical in- 
teractions in star clusters (Sigurdsson & Phinney 1993) and 
thus enhance the total number of BH-BH mergers in the 
universe. There are two types of interaction that create BH- 
BH binaries. The first is few-body binary formation where 
a close encounter between multiple stars transfers kinetic 
energy to one star at the expense of the relative energy be- 
tween two others. This leaves a bound pair and an escaper. 
The second is exchange where a single star interacts with a 
binary and is exchanged into the binary in place of one of 
the original members. Due to equipartition of energy both 
of these processes favour the escape of the lightest member 
of the interaction (Hills & FuUerton 1980; Heggie, Hut & 
McMillan 1996) and thus preferentially introduce massive 
objects, such as BHs, into binary systems. Another impor- 
tant effect is binary hardening where "hard" binaries (bi- 
naries with a binding energy greater than the kinetic ther- 
mal energy of the cluster) can have their binding energy in- 
creased and thus their period shortened (Heggie 1975) . This 
process reduces the separation between the binary members 
and thus can bring a BH-BH binary into the gravitation 
wave regime more quickly than would isolated evolution in 
galactic fields. Interactions between binaries can result in 
both exchange and binary hardening for both binaries. Fi- 
nally interactions can destroy binaries if the kinetic energies 
of the centres of mass are sufficiently high and thus BH-BH 
binaries can also be removed by dynamical interactions. All 
of these processes are enhanced by high stellar densities and 
thus proceed most rapidly in the very concentrated cores 
of star clusters. Equipartition of energy requires that the 
most massive objects in a star cluster sink to the centre in 
a process called mass-segregation (Spitzer 1987). Since BHs 
and binaries rapidly become the most massive objects in the 
system they will rapidly migrate to the centre of the cluster 
where interaction rates are highest. This means that BHs 
will tend to experience a disproportionately large number 
of interactions and this further enhances the effect of star 
cluster dynamics on the BH-BH binary population. 

The effect of star cluster dynamics have been studied 
before using limited dynamical models. Giiltekin, Miller & 
Hamilton (2004) and O'Leary et al. (2006) have investigated 
the BH-BH binary population in star clusters assuming the 
BHs form completely mass-segregated subsystem that inter- 
acts only with itself. They find that the BHs and BH-BH bi- 



naries interact strongly with each other leading to both the 
formation and destruction of BH-BH binaries. Depending 
on the initial conditions assumed for their clusters, O'Leary 
et al. (2006) find ~ 1 — 10 BH-BH merger detections per 
year using the parameters for advanced LIGO. Up to 70% 
of these mergers occur in BH-BH binaries that have been dy- 
namically ejected from the clusters. These simulations con- 
tained no treatment of stellar evolution. By contrast Sad- 
owski et al. (2008) have performed cluster simulations with 
stellar evolution but assuming that the BHs and BH-BH bi- 
naries always remain in dynamical equilibrium with the rest 
of the cluster. This is in some sense the opposite dynami- 
cal assumption made in O'Leary et al. (2006) and leads to 
a much higher merger rate of ~ 25 — 3000 detections per 
year, depending again on the initial conditions. This is be- 
cause although there are fewer interactions that form BH- 
BH binaries, once formed such a binary is unlikely to interact 
with another BH and be disrupted. They claim only ~ 10% 
of mergers occurred outside the clusters. Finally Portegies 
Zwart & McMillan (2000) have performed direct A/'-body 
simulations of star clusters focusing on BHs and the BH-BH 
binary population. These simulations include full cluster dy- 
namics but contain only ~ 4096 particles and as such are 
much smaller than real globular clusters. Portegies Zwart & 
McMillan (2000) find a detection rate of a few detections per 
year to almost one a day depending on the assumed globular 
cluster formation history and the composition of the globu- 
lar cluster population. They also find that the vast majority 
of the BH-BH binaries are ejected from the cluster and thus 
most mergers occur in the field. Consequently their results 
are more consistent with those of O'Leary et al. (2006) than 
Sadowski et al. (2008). 

In Downing et al. (2010) (hereafter Paper I) we have 
revisited this problem using a Monte Carlo star cluster sim- 
ulation code that includes stellar evolution, a fully self- 
consistent treatment of the global cluster dynamics, and 
is capable of simulating a number of particles compara- 
ble to real globular clusters. We demonstrate that the BHs 
strongly mass-segregate and thus confirm that the approxi- 
mation made in O'Leary et al. (2006) is favoured over that 
of Sadowski et al. (2008). We found several potential LISA 
sources in our simulations, four of which would be detectable 
in nearby globular clusters. Our results differ from O'Leary 
et al. (2006) in that we find no BH-BH mergers within our 
clusters, possibly due to our more approximate treatment 
of few-body interactions. This result is, however, consistent 
with the direct A'^-body results of Portegies Zwart & McMil- 
lan (2000). In Paper I we found that many BH-BH binaries 
were ejected from the cores of our clusters with very high 
binding energies, consistent with the results of both Porte- 
gies Zwart & McMillan (2000) and O'Leary et al. (2006). 
These binaries are massive enough and have short enough 
periods that they should merge in a galactic field within 
a Hubble time (Th)- Investigating these escaping binaries 
is the subject of this paper. In § 2 we briefly describe our 
Monte Carlo code and initial conditions. In § 3 we describe 
the properties of the escaping binaries. In § 4 we consider the 
merger rate of escaping binaries and we translate this into a 
detection rate for ground-based gravitational wave detectors 
in § 5. In § 6 we consider whether any of our binaries can 
be detected by the planned LISA space-based detector. We 
discuss our results in § 7 and conclude in § 8. 
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2 METHODS AND SIMULATIONS 

Here we briefly outline our numerical methods and initial 
conditions. We base our results on the simulations described 
in Paper I and the interested reader is referred there for 
further details. 



2.1 The Monte Carlo Code 

We simulate star clusters using a Henon-type Monte Carlo 
code (Henon 1971) with improvements for both global and 
binary dynamics developed by Stodolkiewicz (1982) and 
Stodolkiewicz (1986) as incorporated by Giersz (1998). In 
this code star clusters are assumed to be spherically sym- 
metrical with their global dynamical evolution governed by 
two-body relaxation. Using this approximation the orbits of 
individual stars can be described as plane rosettes defined 
by their energy, E, and angular momentum vector, J. The 
dynamical evolution of the cluster over a time At can then 
be calculated using an appropriate scattering angle chosen 
by Monte Carlo sampling from the theory of two-body re- 
laxation. This method has the advantage of requiring only 
a fixed number of operations per particle and thus scaling 
as 0{N^~^), where A*' is the number of particles in the sim- 
ulation, as opposed to the 0{N'^~'^) of direct A'^-body sim- 
ulations. This allows us to run simulations with 10''' — 10^ 
particles, similar to the number of stars in real globular clus- 
ters, while still being able to perform large parameter space 
studies. 

Strong few-body interactions are not part of the Monte 
Carlo approximation and must be modelled separately by in- 
corporating cross-sections for such events to occur and pre- 
scriptions to determine their outcomes. For three-body inter- 
actions the probability of a close encounter is calculated ac- 
cording to the cross-sections found in Giersz (2001) and the 
outcomes according the the prescriptions of Giersz (1998). 
The cross-sections for binary-single and binary-binary in- 
teractions are calculated following the method of Giersz & 
Spurzem (2003) with the outcome of binary-single interac- 
tions determined by the formulae given in Giersz (1998). The 
outcomes of binary-binary interactions are based on the nu- 
merical few-body scattering experiments of Mikkola (1984) 
as implemented by Stodolkiewicz (1986). The probabilities 
for exchanges during binary-single and binary-binary inter- 
actions are taken from Heggie, Hut & McMillan (1996). 

The code includes treatments for single and binary stel- 
lar evolution using the SSE (Hurley, Pols & Tout 2000) and 
BSE (Hurley, Tout & Pols 2002) packages (Giersz, Heggie & 
Hurley 2008). These packages include prescriptions for stel- 
lar and binary evolution from the zero-age main sequence 
to the stellar remnant at a variety of metallicities and in- 
clude analytic treatments of mass-transfer and perturbed 
evolution in binary systems. For our purposes, the primary 
effect of metallicity is higher-mass BHs at low metallicity 
due to less efficient line driving of stellar winds on the gi- 
ant and asymptotic giant branch (Belczynski, Kalogera & 
Bulik 2002). A study of the resulting BH mass distributions 
is given in Belcyznski et al. (2006). Prescriptions for na- 
tal kicks in supernovae remnants (Lyne & Lorimer 1994) 
are also included. They are drawn from a Maxwellian dis- 
tribution with a peak centred at ~ 190 km s~^ based on 
the proper motion studies of Hansen & Phinney (1997) and 



then reduced in proportion to the mass of material accreted 
during fallback as calculated in Belczynski, Kalogera & Bu- 
lik (2002). An estimate for the gravitational wave inspiral 
timescale based on the quadrupole approximation (Peters 
1964) is also included. Tidal truncation of the cluster is 
treated using the prescriptions of Baumgardt (2001) (Giersz 
2001; Giersz, Heggie & Hurley 2008). 

The code agrees well with direct A''-body simulations 
in the case of equal-masses (Giersz 1998), a mass function 
(Giersz 2001, 2006) and when stellar evolution has been in- 
cluded (Giersz, Heggie & Hurley 2008). Simulations using 
this code have also been shown to re-produce the physical 
parameters of the observed star clusters M67 (Giersz, Heggie 
& Hurley 2008), M4 (Giersz & Heggie 2008) and NGC 6397 
(Giersz & Heggie 2009). Despite these successes the code 
has some limitations. Of particular interest for our study 
is the Spitzer instability (Spitzer 1987) where massive ob- 
jects, such as BHs, can fall out of energy equipartition and 
form a strongly interacting subsystem in the cluster core. In 
such a situation the assumption of weak scattering may no 
longer hold and the Monte Carlo approximation can break 
down. Heggie & Giersz (2009) have compared the behaviour 
of the core region of Monte Carlo and direct A''-body simula- 
tions of NGC 6397 and have found good agreement between 
both sets of simulations in the binary binding energies and 
escape rates. This indicates that core dynamics are being re- 
produced properly. Artificially small bound subsystems can 
also be detected by sudden drops in the central potential of 
the cluster. None of our simulations demonstrate such an ef- 
fect, indicating that there is no major unphysical behaviour 
in the cluster cores. Finally we showed that the binding en- 
ergy of the escapers was in good agreement with that found 
in O'Leary et al. (2006) (and Portegies Zwart & McMillan 
(2000)) indicating that the energetics of binaries in the core 
are being treated accurately. A more serious concern is the 
lack of direct few-body integration in our binary-single and 
binary-binary interactions. This reduces the number of pos- 
sible outcomes for these interactions and eliminates both the 
possibility of mergers during the interaction and the forma- 
tion of hierarchical triples that can pump up the eccentric- 
ity of the inner binary through resonances (Kozai 1962) and 
enhance gravitational wave emission. Eccentricities are also 
not generated self-consistently in interactions but are choose 
randomly. Since most interactions will tend to favour the 
formation of high eccentricity orbits, this means that our 
eccentricities will be under-estimated. The combination of 
these effects means that our merger rates will be lower lim- 
its and could increase with a more accurate treatment of 
few-body interactions. 

2.2 Initial Conditions 

The initial conditions for our simulations were described in 
Paper I and we merely summarise them here. Our results are 
based on a set of 160 simulations, each with N — 5x 10^ par- 
ticles, a Kroupa power-law IMF (Kroupa, Tout & Gilmore 
1993) with a lower slope of ai = 1.3, an upper slope of 
Q; = 2.3, a break mass of 0.5 Mq and masses between 0.1 
Mq and 150 Mq. The models are initialised with Plummer 
profiles with a tidal radius of rt — 150 pc. We have used 
two initial binary fractions (ft), 10% and 50%, with initial 
binary parameters taken from Kroupa (1995). We use two 
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Figure 1. The number of BH-BH escapers per Gyr summed over all ten independent realisations of each set of initial conditions. The 
error bars gives the rms scatter between the independent realisations. 



metallicities, Z = 0.02 (solar) and Z = 0.001 (low). Finally 
we use four different initial concentrations defined by the 
ratio of the tidal radius to the initial half- mass radius (r;j), 
rt/rh G 21, 37, 75, 180, corresponding roughly to initial num- 
ber densities within the half-mass radius of 10^ 10^ 10*, lO'"^ 
pc~^. Smaller half-mass radii reduce the half-mass relax- 
ation time (Spitzer 1987), 



trh = 0.138- 



1)1/201/2 In-yiV 



(1) 



where (m) is the average mass in the system, G is Newton's 
gravitational constant and 7 — 0.02 is an empirically deter- 
mined constant, in the highest concentration clusters and 
thus leads to faster dynamical evolution. Taken together, 
this gives 16 possible combinations of initial conditions. We 
have then performed ten independent realisation of each sim- 
ulation in order to constrain statistical fluctuations for a 
total of 160 individual simulations. We conflrmed that ten 
simulations is sufficient for accurate statistics. The initial 
conditions are summarised in Table 1 of Paper I. The sim- 
ulations have been performed on the HLRS supercomputer 
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Table 1. The number of BH-BH binaries that escape from the 
cluster and the number of BH-BH mergers. The first column iden- 
tifies the initial conditions. The second column gives the total 
number of BH-BH binary escapers summed over all ten inde- 
pendent realisations of each simulation. The third column gives 
the number of BH-BH escapers averaged over all ten indepen- 
dent realisations. The fourth column gives the number of BH-BH 
escapers that merge within a Hubble time summed over all ten 
independent realisations. The fifth column gives the number of 
BH-BH escapers that merge averaged over all ten independent 
realisations. The uncertainty in columns three and five is the rms 
scatter across the ten independent realisations. 
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Figure 2. The total number of BH-BH escapers per Gyr summed 
over all 160 simulations. 



in Stuttgart. Each is run on a single processor and it takes 
between 4 and 24 hours to simulate a cluster for 1 Th- 



3 PROPERTIES OF ESCAPING BINARIES 

In Paper I we discovered no BH-BH mergers within our 
clusters although we found a few potential LISA sources. 
We found, however, that many of our clusters ejected hard 
BH-BH binaries with binding energies > lOOOfcsT where 
fcflT = McoreCcore/2A'^corG IS the thermal energy in the core 
of the cluster, Mcorc is the total mass of the core, A^coro 
is the number of stars in the core and cTcorc is the veloc- 
ity dispersion in the core. This result is consistent with the 
binding energy distribution found by both Portegies Zwart 
& McMillan (2000) and O'Leary et al. (2006). We specu- 
lated that these hard binaries may merge within a Hubble 
time and be sources of gravitational waves for both ground- 
and space-based detectors. 

In Table 1 we present the number of escapers from each 
of our simulation both summed and averaged over all ten in- 
dependent realisations. The trends in Table 1 follow similar 
trends to those found for the formation of BH-BH binaries 
in Table 3 of Paper I, namely more BH-BH escapers at high 
fb, more BH-BH escapers at higher initial concentration and 
more BH-BH escapers at low metallicity. The trend with 
is due both to the larger number of stars and thus the larger 
number of BHs in high-Jt simulations and to the larger num- 
ber of binaries available for BHs to be exchanged into. The 
trend with initial concentration is also clear; at higher den- 
sity there are more interactions per unit time and hence 
more dynamical BH-BH binary formation and ejection over 
1 Th- The trend with metallicity is due both to the larger 
number of BHs and their greater mass at \ow-Z. The mass- 
segregation timescale governing how quickly massive objects 
sink to the cluster centre is given by: 

tms trh (2) 

mi 

where mi > m2 (Watters, Joshi & Rasio 2000; Khalisi, 
Amaro-Seoane & Spurzem 2007) and the more massive 
BHs and BH-BH binaries in low-metallicity clusters mass- 
segregate faster. Thus the massive BHs enter the high den- 
sity cluster cores more quickly than their lower-mass coun- 
terparts and the dynamical evolution is accelerated. Col- 
umn two of Table 1 shows that that there is very little 
simulation-to-simulation scatter in the number of BH-BH 
binaries ejected. 

Figure 1, which gives the number of BH-BH escapers 
per Gyr relative to the start of the simulation from which 
they escaped, confirms this picture. There are more BH-BH 
escapers at high fb, high initial concentration and low metal- 
licity. Figure 1 also shows the effect the initial conditions 
have on the time and rate at which BH-BH binaries escape 
from the simulations. For low initial concentrations and solar 
metallicity the escape of BH-BH binaries commences late in 
the evolution of the cluster and happens at a slow, fairly con- 
stant rate. By contrast at high initial concentrations and low 
metallicity there is a large burst of escapers very early but 
the escape rate drops-off later in the life of the cluster. This 
is due to the shorter relaxation time and mass-segregation 
timescale in these simulations. In the high concentration, 
low metallicity simulations the BHs and binaries mass seg- 
regate and sink to the high density core more swiftly where 
they interact and are ejected rapidly. This process proceeds 
more slowly in the low initial concentration, high metallicity 
clusters and therefore these clusters produce fewer ejections 
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Figure 3. The period distribution of BH-BH escapers binned uniformly in log space and summed over all ten independent realisations 
of each set of initial conditions. Error bars give the rms scatter between the ten realisations. 



within a Hubble time. The initial binary fraction does not af- 
fect the location of peak BH-BH ejection but does affect the 
overall number of BH-BH binaries ejected and clusters with 
a higher initial binary fraction eject more BH-BH binaries 
than their lower /t, counterparts. In figure 2 we present the 
total number of escapers summed over all 160 simulations 
and binned by Gyr assuming all clusters form at the same 
time. It is apparent that the escape rate will be highest early 
in the life of the cluster population where it is dominated 
by low metallicity, high density clusters but that there will 
be an appreciable number of escapers at late time as well. 



Figure 3 gives the period distribution of all escaping bi- 
naries. All periods are short with few more than a year and 
many less than a day. Binaries with periods of less than a 
few days are capable of merging within a Hubble time due to 
gravitational radiation (Peters 1964) and are thus potential 
sources for ground-based detectors. Furthermore, unlike the 
binaries that remain within the cluster, these binaries will 
not be subject to disruption in the further interactions that 
prevent mergers within the cluster as described in Paper 
I. There is little variation in the period distribution with 
cluster parameters other than a tendency towards shorter 
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Figure 4. The eccentricity as a function of period for BH-BH binaries at the time of escape from the cluster. Each panel shows all 
binaries from all all independent realisations of a single set of initial conditions. From top to bottom: = 0.1 and Z = 0.02, /t, = 0.5 
and Z = 0.02, = 0.1 and Z = 0.001, ft, = 0.5 and Z = 0.001. From left to right n/r^ = 21, 37, 75, 180. 



periods in high density clusters. This is because the core ve- 
locity dispersion and thus the core kinetic energy is higher 
in these clusters meaning that a binary must be more en- 
ergetic and thus shorter period in these clusters in order to 
survive and be ejected. 

In figure 4 we present the eccentricity of escaping BH- 
BH binaries as a function of their period at the time of 
escape from the cluster. High eccentricity can substantially 
enhance the power radiated in gravitational waves (Peters & 
Mathews 1963; Peters 1964) and thus high eccentricity bina- 
ries will have shorter merger times than their circular coun- 
terparts. Eccentric binaries also emit gravitational waves in 
higher harmonics of the orbital frequency (Peters & Math- 



ews 1963) and this can move long-period binaries into the 
LISA band. Figure 4 shows that our simulations produce 
BH-BH escapers with a large range of eccentricities at all 
periods indicating a rather high probability of short-period, 
eccentric binaries. These eccentricities must be treated with 
some caution since they are not self-consistently generated 
but rather the product of a random choice at the end of an 
interaction. In reality, however, interactions tend to produce 
high eccentricity binaries so if anything figure 4 underesti- 
mates the average eccentricity to be expected in the BH-BH 
escapers. As we will demonstrate in § 4 many of the binaries 
have a period and eccentricity combination that will lead 
to a merger within a Hubble time and thus we predict that 
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Figure 5. M^hirp for BH-BH binaries at the time of escape from the cluster (crosses) and BH-BH mergers (solid triangles). Each panel 
shows all binaries in all ten independent realisations for a single set of initial conditions. 



star clusters will enrich galactic fields witli BH-BH merger 
candidates. We also note that some of these binaries may be 
in the LISA band with significant eccentricity upon ejection 
and thus we predict that escaping BH-BH binaries from star 
clusters may be eccentric LISA sources. We will investigate 
this possibility further in § 6. 

Finally in figure 5 we show the chirp masses Mchirp of 
all BH-BH binaries at the time of escape as well as for all 
BH-BH mergers (see § 4). The chirp mass, 



Mchirp 



(mim2)^'^^ 
[nil + 7712)^/^ 



(3) 



where mi is the mass of the primary and 7712 is the mass 
of the secondary, is an important quantity in gravitational 
wave studies since the amplitude of a gravitational wave, ho, 
is proportional to Afchirp (Pierro et al. 2001): 

^5/3^2/3^5/3 

^0 « (4) 



where u is the angular frequency of the binary, r is the dis- 
tance from the binary to the observer, and c is the speed of 
light. Thus it is this quantity rather than the total mass that 
is important for gravitational wave detection. Mchirp is larger 
for the low metallicity systems due to the afore mentioned 
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Figure 6. Merger rates for BH-BH escapers binned per Gyr. Each panel gives the merger rate for one choice of initial conditions summed 
over all ten independent realisations. The error bars give the rms scatter between the independent realisations. 



higher mass of BHs derived from lower metallicity progeni- 
tors. It is also clear that in the low metallicity cases the bina- 
ries with the highest chirp mass escape first. This is another 
consequence of mass segregation since the high mass BH-BH 
binaries in the core will interact and be ejected more rapidly 
than their lower-mass counterparts. This is not apparent in 
the high metallicity cases since there the spread in mass is 
not sufficient for this efl'ect to be important (see Paper I). 
It is worth noting, however, that the time-dependence of es- 
caper chirp mass does not translate into a time-dependent 
merger chirp mass. This is due to the time lag between the 



escape of a BH-BH binary and its merger in the field since 
a short period, low mass binary may still merge before a 
long-period, high mass binary. This will be explored in § 4. 



4 BH-BH MERGERS 

Here we calculate the number of escaping BH-BH binaries 
that will merge within a Hubble time and when the merger 
will occur. We use the formalism presented in Peters & 
Mathews (1963) and Peters (1964) in order to calculate the 
gravitational wave inspiral timescale for BH-BH binaries in 
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Time [Gyr] 



Figure 7. The number of BH-BH merger binned per Gyr and 
summed over 160 all realisations of all initial conditions. 



a galactic field. In this approximation the field equations are 
linearised and the gravitational wave radiation is then cal- 
culated in a multipole expansion. The first non-zero order in 
this expansion is the quadrupole term, which can be used to 
calculate the energy and angular momentum carried away 
from the binary by gravitational waves. This can then be 
used to calculate the effect of gravitational radiation on the 
orbital elements of the binary. The orbit-averaged change 
in semi-major axis, a, and eccentricity, e, in the quadrupole 
approximation are: 



(d) = 
and 



64 mim2(mi -I- 7712) 
a3(i_e2)7/2 



304 mim2{mi + m2) 



^ ^ 73 2 ^ 37 4 

1 H e H e 

24 96 



V 304 ; 



(5) 



(6) 



(Peters & Mathews 1963; Peters 1964). For each escaping 
binary we solve equations 5 and 6 using a 4*-order Runge- 
Kutta integrator with time steps chosen such that there 
is never more than a 1% change in a. We integrate from 
the time of escape from the cluster until either the binaries 
merge or 1 Th is reached. This gives us both an inspiral 
timescale for the binary and, after each timestep, a self- 
consistently generated, quasi-static set of orbital parameters 
that can be used to calculate a gravitational wave signal for 
LISA using the formalism of Pierro et al. (2001) as will be 
discussed in more detail in § 6. 

In table 1 we show the number of escaping BH-BH bi- 
naries that merge within a Hubble time. The total number 
of BH-BH mergers tracks the total number of BH-BH es- 
capers however the proportion of escaping BH-BH binaries 
that merge within a Hubble time increases in clusters with 
low metallicity and high initial concentration. The correla- 
tion with metallicity is again due to the larger mass of BHs 



in the low metallicity simulations and the correspondingly 
stronger gravitational wave radiation. The correlation with 
initial concentration is due to the fact that binaries tend 
to have shorter periods in clusters with a high concentra- 
tion and thus have a shorter distance over which to inspi- 
ral and merge. Again there is relatively little simulation-to- 
simulation scatter in the number of mergers. 

Figure 6 shows the number of mergers for each set of 
initial conditions per Gyr where the time of merger of a 
BH-BH binary is given relative to the start of the simula- 
tion that produced it. In this plot we present a sum rather 
than an average over each set of initial conditions. This is 
partly due to the small number of mergers but also because 
once the BH-BH binaries escape from their parent clusters 
they will form a single population in a galactic field. The 
overall shapes broadly follow those found for the number of 
escapers per Gyr in figure 1 but with a time delay between 
when escapers and mergers begin and when the peak number 
of escapers and mergers per Gyr is found. The mergers are 
delayed because the inspiral in a galactic field is not an in- 
stantaneous process and it often takes a binary several Gyrs 
to merge. Again assuming that all star clusters form at the 
same time early in the life of the universe then the merger 
rate in the early universe will be dominated by BH-BH bina- 
ries originating from dense, low-metallicity clusters whereas 
at the current age of the universe clusters of all types will 
be contributing mergers in similar numbers. In figure 7 we 
present the total number of mergers per Gyr summed over 
all realisations of all initial conditions. The profile is flat- 
ter than for the total number of escapers in figure 2 but still 
shows a peak in the merger rate while the cluster population 
is young. This implies that if we can detect BH-BH mergers 
at large redshift then mergers from young clusters in the 
early universe will dominate the detection rate. It also indi- 
cates that if dense, low-metallicity clusters form later in the 
universe they will be major contributors to the local BH-BH 
merger rate. 

Figure 8 give the inspiral timescale for each BH-BH 
binary from the time of escape from the cluster until the 
merger as calculated by equations 5 and 6 and as a function 
of time of escape from the cluster. There are no particular 
trends in inspiral time with escape time in any of the simula- 
tions. This implies that there is no clear link between merger 
time and escape time and that it is probably impossible to 
derive an escape rate from an observed merger rate. It also 
explains why there is no evolution in Mchirp of mergers as a 
function of time. High Mchirp binaries with a relatively long 
period that escape early are as likely as low A/chirp binaries 
with a short period that escape late. Thus the trend in the 
early escape of high Mchirp binaries does not translate into 
a similar trend in the merger rate. 



5 DETECTION RATES FOR GROUND-BASED 
DETECTORS 

Now we must determine if the BH-BH mergers that are pro- 
duced by our simulations can be detected by gravitational 
wave observatories. We follow the approach of O'Leary et al. 
(2006) who estimate the detectability of a BH-BH merger by 
comparing it to a template merger with a known signal-to- 
noise ratio at a given distance. A merger with a red-shifted 
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Figure 8. The merger timescale as a function of escape time for all escaping BH-BH binaries that merge. Each panel gives merger 
timescales for all ten independent realisations of each set of initial conditions. 



chirp mass of TVfchirp = (1 + 2m)Mchirp, where Zm is the red- 
shift of the merger, is considered detectable at a luminosity 
distance Dl = (1 + z,n)Dprop, where Dprop is the proper 
distance to the merger, if (O'Leary et al. 2006): 

Dl,0 ( McMrp Y''' 1 s(.foft) , . 

Dl ^Mchi.p,o J V s(/off.o) 

where Dl,o is the luminosity distance at which a merger of 
red-shifted chirp mass A1chirp,o is known to be detectable 
at a given signal to noise ratio. For the planned advanced 
LIGO detector it is estimated that a equal-mass, NS-NS 



binary with il'/chirp = 1.2Mq can be detected with a signal- 
to-noise ratio of 8 at a distance of 191 Mpc (O'Shaughnessy 
et al. 2005; Harry et al. 2005). The distance for the current 
generation of detectors is some ten times lower. The detector 
response function, s(/off) is given by (O'Leary et al. 2006); 

rfa!! (rl\-7/S 

sif..)^l ^-ll^df (8) 

where Sjv(.f') is the noise spectrum of the detector and is 
approximated as (Cutler & Flanagan 1994; O'Leary et al. 
2006): 
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Figure 9. Detection rates by cluster type. Each panel gives the expected detection rate if the entire cluster population in the universe 
was composed of identical clusters, each with the corresponding initial conditions. The x-axis gives the lookback time to Tform in Gyr. 
The solid line with crosses is for Dhfi = 19.1 Mpc, the dotted line with triangles is for Dl.o = 191.0 Mpc and the dashed line with 
squares is for -Dj^^o = 1910.0 Mpc. 



Siv(/') cx 
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(9) 



with /o = 70 Hz. A constant of proportionality is not re- 
quired in equation 9 because it will cancel out in equation 7. 

The cut-off frequency, /off is the frequency at which the 
merger occurs and detailed relativistic modelling becomes 
necessary. For the reference inspiral we use the estimate of 
Cutler & Flanagan (1994) which yields /off,o ~ 720 Hz for 



the NS-NS reference binary we have chosen. For the BH-BH 
inspiral the frequency is generally taken to be that of the 
last circular orbit before the final plunge. We use the same 
circular estimate as O'Leary et al. (2006): 



/off 



200 

V M 



Hz 



(10) 



where M is the total mass of the inspiralling binary. 

From our simulations we have Mchirp for each of our 
mergers and the time of the merger relative to the start of 
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the simulation. By assuming a formation time, Tform, for 
each simulation relative to the age of the universe we can 
calculate a lookback time, redshift and proper distance to 
the merger and use equation 7 to determine if the merger can 
be detected. To convert this to a merger rate per year for 
ground-based detectors we must make assumptions about 
the formation history and density of globular clusters in the 
universe. For simplicity we assume that all globular clusters 
form in a single burst at some time Tform in the past. Again 
for simplicity we assume a constant number density of no ~ 
SAh^ clusters Mpc"'^. This estimate comes from Portegies 
Zwart & McMillan (2000) and is based on observations of 
the local universe taking into account the different specific 
frequencies of clusters in galaxies of different Hubble types. 
We then bin the detections in uniform bins in lookback time, 
calculate the proper distance, -Dprop.i, to the edge of each 
bin, and use the equation (O'Leary et al. 2006) 



M-^ N^4n 3 



d: 



prop,l~ 1 



)(l + 20" 



(11) 



to calculate the detection rate. Ni is the number of detec- 
tions calculated from our simulation in bin i, At is the bin 
width and the factor {1 + Zi)~^ comes from the cosmological 
time dilation of the detection rate. We use the cosmologi- 
cal parameters h = 0.72, flm = 0.27 and flLambda = 0.73. 
We note that it would be possible in principle to perform a 
more careful estimate of the density of globular clusters in 
the universe by using press-schechter formalism and making 
assumptions about the baryon fraction per halo, the fraction 
of baryons that become globular clusters, the evolution of 
the globular cluster mass function and globular cluster for- 
mation rate. This would, however, introduce a large number 
of uncontrolled parameters into our modelling and since our 
simulations do not cover the whole range of globular cluster 
masses convolving them with the evolution of the globular 
cluster mass function seems unjustified. For these reasons we 
will limit ourselves to simple assumptions about the globu- 
lar cluster population in the universe and save more detailed 
modelling of the cluster population for when we have a larger 
set of simulations. 

In figure 9 we present the detection rate per year as a 
function of lookback time to the assumed formation event. 
We have chosen a lookback time to Tform to be between 1 and 
13 Gyr. We have also calculated rates for three different de- 
tector sensitivities (parametrised by Dl,o)- The dotted line 
represents advanced LIGO with the solid line representing 
a detector ten times less sensitive (corresponding roughly 
to the detectors currently in operation) and the dashed line 
corresponding to a detector ten times more sensitive. For the 
current generation of detectors there is a very low probabil- 
ity of detection since the rate seldom exceeds 0.1 detection 
per year and is often much lower. The current detectors are 
not sensitive enough to detect mergers out to cosmological 
distances and thus probe a volume that is too small to pro- 
duce a substantial detection rate. By contrast the prospects 
for advanced LIGO are much better with most of the high- 
density, low metallicity clusters promising detection rates of 
more than 1 yr^^ and some of the most optimistic provid- 
ing ~ 50 detections yr~^. The advanced LIGO detection rate 
demonstrates some interesting trends with formation time. 
The low density, high metallicity simulations only produce 




Figure 10. The detection rate assuming all types of clusters are 
present equally in the universe with an overall number density of 
no = 8Ah^ Mpc~^. The solid line with crosses is for Dj^ q = 19.1 
Mpc, the dotted line with triangles for Dj^ q = 191.0 Mpc and 
the dashed line with squares is for D^ o = 1910.0 Mpc. 



detections if the clusters are assumed to form very early 
in the universe. This is because these clusters evolve slowly 
and do not produce BH-BH mergers until late in their lives 
(see figure 6). By contrast the higher density, lower metal- 
licity clusters can provide a significant detection rate even 
if they form quite late in the universe. Indeed these clusters 
produce a peak detection rate for a Tform of only 5-6 Gyr 
ago since, according to figure 6, they evolve more quickly 
and produce the majority of their mergers while still young. 
Thus they have the potential to contribute significantly to 
the merger rate while only a few Gyrs old. The peak at 5-6 
Gyr comes form the balance between detection volume and 
detectability. For formation times less than 5-6 Gyr ago ad- 
vanced LIGO will be able to detect all mergers and thus by 
choosing an earlier formation time the volume of space in 
which globular clusters exist and thus the number of glob- 
ular cluster increases, and thus the detection rate goes up. 
Clusters that form more than ~ 6 Gyr ago are sufficiently 
distant that advanced LIGO cannot detect all of the merg- 
ers and in particular the early peak merger rate moves out 
of range for LIGO detection. Thus the detection rate starts 
to fall again. For a ten times greater sensitivity (currently 
well beyond our technical capabilities) we would be able to 
detect all mergers at much larger distances and thus the de- 
tection rate continues to increase with detection volume up 
to a Tform of ~ 10 Gyr ago. 

Figure 10 gives the detection rate assuming that all of 
our simulated clusters are present in the universe in equal 
numbers. The shape is very similar to that of the low metal- 
licity, high density clusters, indicating that it is these clus- 
ters that will dominate the detection rate. The overall detec- 
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tion rate is, however, much lower. The detection rate for the 
current generation of detectors drops to little more than 0.01 
per year and for advanced LIGO the rate drops to a peak of 
~ 4 per year. This is because the clusters that produce few 
mergers serve only to dilute the detection rate. 4 BH-BH 
detections per year is, however, still twice that predicted by 
Belczynski et al. (2007) for galactic field populations and 
thus star clusters still provide a modest enhancement to the 
BH-BH detection rate. 



6 DETECTION RATES FOR SPACE-BASED 
DETECTORS 

Most of the escaper binaries that are potential LISA sources 
will be in the LISA band for period of time that is less than 
the spread in ages of Galactic globular clusters. However, 
if a source is in the LISA band, it will appear as a quasi- 
monochromatic source for the duration of the LISA mis- 
sion. Therefore, we will express the detectability of escaper 
BH binaries in terms of the probability that a source will 
be in band during the LISA mission. We do this by dis- 
tributing the population of escaper BH binaries throughout 
the Galactic halo and determining the lifetime of each de- 
tectable source. Using the spread in globular cluster ages, 
we can then express the probability of detection as a ratio 
of the lifetime to the spread in globular cluster ages. For this 
paper, we adopt a spread of 4 Gyr. 

In order for a binary to be a LISA source, it must lie 
within the LISA sensitivity band during the expected age 
range of the Galactic globular clusters. We take the age 
range to be between 10 and 14 Gyr (Forbes & Bridges 2010). 
We determine the if a binary is potentially within the LISA 
band by considering whether the harmonic with maximum 
power has a frequency above the bottom of LISA's sensitiv- 
ity band, which we take to be 5 x 10~^ Hz. Out of all 160 
simulations, there are 91 escaped binaries within the LISA 
band during the 10-14 Gyr time- frame. The number of es- 
capers within the LISA band are binned by Gyr and shown 
in figure 11. 

In order to distribute the escapers throughout the 
Galactic halo, we need to adopt a density distribution for 
the escaped binaries. Because we have not explicitly mod- 
eled individual globular clusters, we do not trace the path of 
the escaped binaries through the Galactic potential. Instead, 
we note that the average binding energy of the escapers can 
give an estimate of their recoil velocity following the last 
encounter. For an upper bound on the recoil velocity, we as- 
sume that the binding energy of the escaped binary is equal 
to its recoil kinetic energy, while for a lower bound, we follow 
Portegies Zwart & McMillan (2000) and assume that the re- 
coil energy is 1/15 of the final binding energy. We show the 
number binaries binned by 20 km/s intervals for all escapers 
within the LISA band in figure 12, using both prescriptions 
for the recoil energy. We note also that the escape veloc- 
ity distribution is fairly fiat over the lifetime of the cluster, 
with a slight trend toward lower velocities at later times (al- 
though the spread in velocities increases). This is shown in 
figure 13 for both the upper and lower bounds of velocities. 

The lower bound recoil velocities are comparable to the 
escape velocity for a typical globular cluster, and conse- 
quently the escapers are most likely to be found at Galacto- 




Time [Gyr] 

Figure 11. The total number of escape times for LISA source 
BH-BH escapers per Gyr summed over all 160 simulations. 




320 



Figure 12. The total number of recoil velocities for LISA source 
BH-BH escapers per 20 km/s interval. The solid line is the upper 
bound, assuming 50% of the initial binding energy goes into recoil 
energy during the last interaction. The lower bound, using the 
prescription of Portegies Zwart and McMillan Portegies Zwart & 
McMillan (2000) is the dashed Une. 



centric radial distances that are comparable to the radial 
distribution of the globular cluster system. We use the Har- 
ris catalog Harris (1996) to obtain these distances. On the 
other hand, the upper bound recoil velocities are compara- 
ble to the Galactic rotational velocities at globular cluster 
distances. In this case, the escapers are likely to be found 
distributed according to the Galactic halo distribution. We 
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Figure 13. The estimated escape velocity plotted against the 
time of escape for both upper and lower bounds. 



use a simplified spherical density profile (Zinn (1985); Mor- 
rison (1996); Siegel et al. (2002)): 

Phalo oc f 1 + — ) (12) 
\ ao / 

where ao = 3.5 kpc is the scale radius of the halo. We impose 
a cut-off radius of Vco = 300 kpc 

Using a Monte Carlo technique, we place each poten- 
tially detectable BH-BH binary in 100 locations determined 
by the globular cluster distribution for an upper bound on 
the detection likelihood and 100 locations determined by 
the halo distribution for the lower bound. The initial phase 
and orientation are randomly chosen at each location. For 
each location, we compute the orbital period and eccentric- 
ity when the binary first becomes detectable in the LISA 
data stream. We determine detectability by imposing a de- 
tection threshold on the signal-to- noise ratio (p), computed 
according to: 



\Kf)\' 

Snif) 



■df 



(13) 



where h{f) is the Fourier transform of the response of LISA 
to the gravitational wave and S„(/) is the power spectral 
density of the expected noise in LISA. We include both in- 
strument noise and an estimate of the Galactic white dwarf 
binary foreground from Ruiter et al. (2010). We compute the 
response of LISA using the time domain as described in Be- 
nacquista, DeGoes & Lunder (2004), taking the barycen- 
tered waveform for eccentric binaries from Pierro et al. 
(2001), carried out to the n ~ 1300 harmonic. We set a de- 
tection threshold of p ^ 10 in a single interferometer (which 
corresponds to a combined p ^ 14 in two channels of the 
LISA data stream). 

Once the initial orbital period and eccentricity at the 
onset of detection have been determined, we compute the 
lifetime of the binary in the LISA band using Peters (1964) 
to find the time to coalescence. We note that the actual 
coalescence will occur above of the LISA band, but that 



Table 2. The average lifetime and likelihood for LISA detection 
for escapers using Globular cluster radial distribution or halo ra- 
dial distribution. Errors are one standard deviation obtained from 
100 realisations of each distribution. 
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50sol37 
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50sol75 
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31 ± 
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33 ± 


36 
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19 ± 


51 
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257 ± 230 


64 ± 


57 


146 ± 331 


37 ± 


83 
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93 ± 133 


23 ± 


33 


51 ± 171 


13 ± 


43 


lOlowlSO 


187 ± 154 


47 ± 


39 


108 ± 226 


27 ± 


57 


501ow21 


40 ± 77 


10 ± 


19 


24 ± 121 


6 ± 


30 


501ow37 


400 ± 308 


100 ± 


77 


225 ± 420 


56 ± 105 


501ow75 


333 ± 214 


83 ± 


54 


189 ± 284 


47 ± 


71 


501owl80 


405 ± 234 


101 ± 


59 


227 ± 319 


57 ± 


80 



the time spent in this phase is negligible compared with the 
time spent in the LISA band. Thus, for each BH-BH binary, 
we have 100 possible lifetimes in the LISA band, which we 
condense into a mean lifetime and its associated variance. In 
order to estimate the likelihood of a given model producing 
a detectable BH-BH binary from its population of escapers, 
we average over the lifetimes of each BH-BH binary in each 
model. We also average over the variances to obtain a spread 
in the lifetimes. We do this for both the globular cluster 
distribution and the halo distribution. The results are shown 
in Table 2. We also show likelihoods in Figure 15. We note 
that the models with low metallicity and high initial binary 
fraction have the highest likelihood of ejecting a BH-BH 
binary that will be detectable with LISA. However, even in 
the most optimistic case with every Galactic globular cluster 
represented by the highest likelihood model (501ow37), the 
likelihood of LISA detection is less than 3 %. Thus, detection 
of a BH-BH binary produced in a Milky Way globular cluster 
is unlikely. 

The majority of BH-BH binaries in the LISA band are 
close to circular with very little model-to-model variation. 
There are, however, a few binaries with very high eccentric- 
ity and, as is shown in Figure 14, these are the ones that 
have recently escaped from their parent clusters. This result 
is not surprising since the binaries that have been in the field 
for the shortest period of time are closer to their last inter- 
action and thus to the last chance for signficant eccentricity 
to be imparted before gravitational wave radiation carries it 
away. This also agrees with the results of Paper I where all 
the LISA binaries that remain in the clusters and are subject 
to frequent interactions remain highly eccentric. This opens 
a small (but non-zero) possibility of a highly eccentric stellar 
mass detection in the LISA band. 

Since BH-BH binaries are strong sources of gravita- 
tional radiation, they can be observed at distances substan- 
tially greater than the scale of the Milky Way. We can esti- 
mate the maximum distance to which a BH-BH binary can 
be detected after an observation time T using: 



(14) 
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Figure 14. The eccentricity of BH-BH binaries in the LISA band as a function of the time they have been in the galactic field for both 
the globular cluster and the halo distributions. 




Figure 15. Likelihood of LISA detection of an escaped BH-BH pj^^^^ Cumulative histogram showing the fraction of bina- 

binary from a single globular cluster for each model. Model num- ^j^^ ^j^^ LISA sensitivity band that arc detectable at a distance 

bers correspond to the order in which models are listed in Table 2. 



where he is the characteristic strain at a given gravitational 
wave frequency, p is the threshold snr for detection, and \/iS'ti 
is the strain spectral density at that frequency. For h^, we 
use (Benacquista 1999): 

, r;5/3 , . 5/3 r2/3 

From the LISA Sensitivity Curve Generator available at 
http : //www. srl . caltech. edu/~shane/ sensitivity/ (Lar- 
son, Hellings, & Hiscock 2002) , we find the minimum strain 
spectral density to be y/S^ — 1.2 x 10^^'Vv^iz at 
/gw = 7 mHz. Setting a detection threshold of p = 10, we 
find that the median maximum distance for all detectable 
binaries to be rmax ~ 7.5 Mpc. The cumulative histogram 



of Tmax for all detectable binaries in all models is shown in 
Figure 16. 



7 DISCUSSION 

There are as yet no detections of gravitational waves in 
either the high-frequency ground-based or low-frequency 
space-based regimes. We can, however, compare our ground- 
based results to previous studies of BH-BH binaries in glob- 
ular clusters using different dynamical assumptions and to 
some tentative predictions for stellar-mass LISA sources in 
the Galactic field. 

As previously mentioned, the investigations of O'Leary 
et al. (2006) and Sadowski et al. (2008) are some of the 
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most complete available and make opposing dynamical as- 
sumptions regarding mass-segregation of BHs and thus the 
interaction rate in the core. In the work of O'Leary et al. 
(2006) BHs are assumed to mass segregate completely and 
form a strongly interacting sub-system in the core of the 
cluster. This leads to a large number of interactions that 
can both create and destroy BH-BH binaries. By contrast 
Sadowski et al. (2008) assume that the BHs do not mass- 
segregate and remain in thermal equilibrium with the rest 
of the stars in the cluster. This leads to a much lower in- 
teraction rate between the BHs and, consequently, a lower 
rate of both BH-BH formation and destruction. In practise 
the lower destruction rate wins and the simulations of Sad- 
owski et al. (2008) predict a yearly detection rate at least 
an order of magnitude higher than that of O'Leary et al. 
(2006). This implies that the predictions of O'Leary et al. 
(2006) represent a lower limit on the BH-BH merger de- 
tection rate while those of Sadowski et al. (2008) represent 
an upper limit. O'Leary et al. (2006) also predict that up 
to 70% of BH-BH mergers will occur in binaries that have 
been ejected from the cluster whereas Sadowski et al. (2008) 
predict that ~ 90% of mergers occur in binaries remaining 
within the clusters. 

The predictions of our simulations, using a fully self- 
consistent treatment of stellar dynamics, agree with those 
of O'Leary et al. (2006) better than those of Sadowski et 
al. (2008). We confirm that few, if any, mergers should take 
place within the clusters and the peak detection rates pre- 
dicted for our denser clusters in Figure 9 fall in the same 
range as those given in Table 2 of O'Leary et al. (2006). 
This agreement is in spite of the fact that both O'Leary 
et al. (2006) and Sadowski et al. (2008) use a similar and 
much more detailed treatment of few-body interactions than 
is present in our Monte Carlo code. For these reasons we 
favour the approximations made by O'Leary et al. (2006) 
rather than those made by Sadowski et al. (2008) and pos- 
tulate that a proper treatment of global globular cluster dy- 
namics may be more important in predicting event rates 
for ground-based detectors than the detailed microphysics 
of the interactions. 

The most significant difference between our predictions 
and those of O'Leary et al. (2006) is that we do not find 
any mergers between BH-BH binaries that remain in our 
clusters, only in the escapers. As we demonstrated in Pa- 
per I, this is because in our clusters BH-BH binaries tend 
to either be disrupted or ejected by dynamical interactions 
before they have a chance to merge due to gravitation wave 
radiation. It is possible that this effect is real, as suggested 
by the results of Portegies Zwart & McMillan (2000), how- 
ever it may also be due to our more approximate treatment 
of few-body interactions. As discussed in Paper I we may be 
overestimating both the disruption rate in binary-binary in- 
teractions and the velocity kick applied to the centre-of-mass 
of each member of an interaction while at the same time un- 
derestimating the eccentricities of binaries formed through 
these interactions. This combination of effects would lead 
to shorter disruption timescales, larger escape probabilities 
and longer gravitational wave inspiral timescales than would 
be expected in a real cluster. Thus it is almost certain that 
our merger rates, particularly those within the cluster, are 
lower limits. A version of our Monte Carlo code with a more 
detailed treatment of few-body interactions is currently in 



development and will allow us to investigate our BH-BH 
escape rates in more detail. 

We have no direct comparisons in the LISA band since 
all studies of globular clusters as sources of resolved stellar- 
mass LISA sources have focused only on those binaries that 
are retained by the cluster (e.g. Benacquista (2001a); Benac- 
quista, Portegies Zwart, & Rasio (2001b); Kocis, Caspar, & 
Marka (2006); Ivanova et al. (2006); Willems et al. (2007)). 

Much more attention has been paid to field binaries. 
Both Nelemans, Yungelson & Portegies Zwart (2001) and 
Belczynski, Benacquista & Bulik (2008) have investigated 
the possibility that BH-BH binaries in the Calactic field 
could be resolved stellar-mass LISA sources. Belczynski, Be- 
nacquista & Bulik (2008) further compared two different 
models for BH-BH binary formation in the Calactic field. 
They found that if binaries can survive the Hertzsprung gap 
without merging then at the present time there may be ~ 10 
resolvable double compact binaries in the Calactic field, half 
of which may be BH-BH. By contrast if most binaries merge 
during the Hertzsprung gap then there would only be ~ 2 re- 
solvable double compact binaries in the Calactic field, none 
of which are expected to be BH-BH. Taking our estimate 
that our 4 Cyr "snapshot" represents the current Calactic 
globular cluster system at face value and noting 85 out of 
our 91 LISA binaries are resolvable we find ~ 1 x 10~^ re- 
solvable BH-BH LISA sources per initial 10'' M0 in globular 
clusters as opposed to the ~ — 1 x 10~^ per initial lO^M© 
in Calactic stellar mass found by Belczynski, Benacquista & 
Bulik (2008). Thus globular clusters do seem to efficiently 
produce LISA sources. We note, however, that it is very un- 
likely that all 85 of the resolvable BH-BH binaries we find 
will be "on" at the same time, hence our low detection prob- 
ability compared to Belczynski, Benacquista & Bulik (2008) 
who's BH-BH binaries will all be on at the same time by con- 
struction. Belczynski, Benacquista & Bulik (2008) concluded 
that there is a reasonable probability that LISA will resolve 
Calactic double compact binaries and that the presence of 
BH-BH detections can be used to constrain the behaviour of 
binaries during the Hertzsprung gap. Our results do not con- 
tradict the first point, although Calactic globular clusters 
may not contribute much to the LISA detection rate. They 
do, however, imply that even if a BH-BH binary is detected 
in the LISA band it is not necessarily evidence that binaries 
survive the Hertzsprung gap phase without merging since 
the binary could be produced by dynamical interactions in 
a star cluster. Our results also indicate that if a highly ec- 
centric BH-BH binary is detected by LISA is likely to have 
experienced a dynamcial interaction in its recent past and 
thus is likely to have originated in a star cluster. Note that 
the converse is not true, a circluar BH-BH binary in the 
LISA band may be either a cluster binary that has existed 
in the field for long enough to cirularise or it may have a 
galactic origin. 

The discovery that star clusters produce BH-BH bina- 
ries that are visible to LISA up to Mpc distances is very 
interesting since it indicates that we may be able to detect 
such binaries in galaxies other than our own. For this rea- 
son we plan to extend our work to other galaxies in our local 
group and determine if this will have a significant effect on 
the probability of LISA detecting a stellar-mass BH-BH bi- 
nary. 



18 J. M. B. Downing, M. J. Benacquista, M. Giersz, & R. Spurzem 



8 CONCLUSIONS 

We have studied the behaviour of the BH population in 
globular clusters, focusing on BH-BH binaries formed by 
dynamical interactions but then inspiralling and merging in 
a galactic field. We find that globular clusters produce such 
binaries quite efficiently and should enhance the BH-BH bi- 
nary detection rate for the next generation of ground-based 
detectors by a modest factor over the rate expected from 
isolated stellar evolution in galactic fields. This enhance- 
ment is in line with the predictions of O'Leary et al. (2006) 
rather than the more extreme factors predicted by Sadowski 
et al. (2008). We find that globular clusters will also pro- 
duce stellar-mass LISA sources resolvable at distances of 
several Mpc. The probability of such a binary being "on" in 
our Galaxy during the three year LISA mission is, however, 
rather low. In future we plan to include few-body integration 
in our code and see if we get better agreement with O'Leary 
et al. (2006) on the number of mergers within clusters. Given 
the large distance to which LISA can detect some of our 
binaries, we are also interested in extending our study to 
other local group galaxies with simulations chosen to more 
closely match the real globular cluster population. We con- 
clude that star clusters can efficiently enrich galactic fields 
with BH-BH binaries and must be taken into account when 
estimating detection rates and population characteristics for 
ground- and space-based gravitational wave detectors. 
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